
clear
clc

cd 'E:\data_replication'

i_k = 2000
KK = 4092

load('estimation/5_supply_side/main_data_supply.mat');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Import Data and Estimates for Product k %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

k = index_price_coeff(i_k,1);    


data_k = data_summary{k};
numProdsTotal = size(data_k,1);
X_MPEC_k = X_MPEC_summary{k};
X_MPEC_opt = X_MPEC_opt_summary{k};
share_k = data_k(:,5);
price_k = data_k(:,7) + 1;
gdp_capita_k = data_k(:,13);
distance_k = data_k(:,14);
population_k = data_k(:,16);
population_log_k = log(population_k);
Y_k = Y_summary{k};
FE_k = FE_summary{k};
market_id_k = data_k(:,end);
T = market_id_k(end,1);
nn = size(Y_k,2);
alpha_constant = X_MPEC_opt(1,1);
alpha_distance = X_MPEC_opt(2,1);
alpha_population = X_MPEC_opt(3,1);
alpha_price = X_MPEC_opt(4,1);
alpha_1 = X_MPEC_opt(12,1);
alpha_FE = X_MPEC_opt((3*(K_A+1) + 1):(3*(K_A+1) + (K - K_A)), 1);

xi_k = X_MPEC_opt((3*(K_A+1) + (K - K_A) + 2 + 1):(3*(K_A+1) + (K - K_A) + 2 + numProdsTotal), 1) + 10;     % Ensure positive values
xi_k = xi_k - alpha_constant - log(price_k)*alpha_price - alpha_population*population_log_k - alpha_distance*distance_k;
q_k = alpha_constant + xi_k + FE_k*alpha_FE;


prodsMarket = zeros(T,1);
prodsMarket_temp = data_k(:,17);
for t = 1:T
prodsMarket_temp2 = prodsMarket_temp(market_id_k==t);
prodsMarket(t,1) = prodsMarket_temp2(1,1);
end

marketStarts = zeros(T,1);
marketEnds = zeros(T,1);
marketStarts(1) = 1;
marketEnds(1) = prodsMarket(1);
for t=2:T
    marketStarts(t) = marketEnds(t-1) + 1;
    marketEnds(t) = marketStarts(t) + prodsMarket(t) - 1;
end

marketForProducts = zeros(numProdsTotal,1);
for t=1:T
marketForProducts(marketStarts(t):marketEnds(t)) = t;
end
numProdsTotal = size(data_k, 1);


delta_k =  alpha_constant + xi_k + log(price_k)*alpha_price + alpha_population*population_log_k + alpha_distance*distance_k;

for t = 1:T   
Y_market = Y_k(t,:);
expmu(marketStarts(t):marketEnds(t),:) = exp(log(price_k(marketStarts(t):marketEnds(t),:))*alpha_1*Y_market + repmat(FE_k(marketStarts(t):marketEnds(t),:)*alpha_FE,1,nn));
end
expmeanval = exp(delta_k);
oo = ones(1,nn);                  
sharesum = sparse(zeros(T,numProdsTotal));  % used to create denominators in logit predicted shares (i.e. sums numerators)
for t = 1:T
    sharesum(t,marketStarts(t):marketEnds(t)) = 1;
end
numer = (expmeanval*oo ).*expmu;        
sum1 = sharesum*numer;
sum11 = 1./(0+sum1);                                                       % Will not result in share_k = EstShare_true
denom1 = sum11(marketForProducts,:);    
simShare_true = numer.*denom1;              
EstShare_true = mean(simShare_true,2);  



% Estimation
%==========================================================================

load('estimation/5_supply_side/companies_k.mat');
companies_k = company_summary{k};

load('data/firm_size_distribution/prod_summary.mat');
prod_pct = prod_pct_summary{i_k};

cd 'E:\data_replication\robustness\firm_heterogeneity'

data_k_summary = [];
f_summary = [];
price_k_summary = [];
mc_summary = [];
q_k_summary = [];
share_k_summary = [];
distance_k_summary = [];
population_log_k_summary = [];
FE_k_summary = [];
m1_summary = [];
N_true_summary = [];
market_id_summary = [];
p_flag_summary = [];

for ttt = 1:T

ttt    

index_market = market_id_k == ttt;
data_market = data_k(index_market,:);
pos_dom = find(data_market(:, 3) == data_market(:, 4));
numProdsTotal_market = size(data_market, 1);
prodsMarket_market = prodsMarket(ttt,1);

marketStarts_market = 1;
marketEnds_market = prodsMarket_market;
marketForProducts_market = ones(prodsMarket_market, 1);

N = table2array(companies_k(:, 4));
N(N < 1) = 1;
N(N > 30) = 30;
N = N + 1;
N = N(index_market);
N_d = N(1,1);
N_d = round(N_d);
N_d_market = ones(numProdsTotal_market, 1);
N_d_market(pos_dom, :) = N_d;

% Extract market-specific data
share_k_market = share_k(index_market,:);
price_k_market = price_k(index_market,:);
q_k_market = q_k(index_market,:);
q_k_market = q_k_market - log(N_d_market);
FE_k_market = FE_k(index_market,:);
population_log_k_market = population_log_k(index_market,:);
distance_k_market = distance_k(index_market,:);
xi_k_market = xi_k(index_market,:);
xi_k_market = xi_k_market - log(N_d_market);
delta_market =  alpha_constant + xi_k_market + log(price_k_market)*alpha_price + alpha_population*population_log_k_market + alpha_distance*distance_k_market + log(N_d_market);
Y_market = Y_k(ttt,:);
expmu_market = exp(log(price_k_market)*alpha_1*Y_market + repmat(FE_k_market*alpha_FE,1,nn));
expmeanval_market = exp(delta_market);
oo = ones(1,nn);                  
sharesum_market = sparse(zeros(1,size(data_market,1)));
sharesum_market(1,:) = 1;
numer = (expmeanval_market*oo ).*expmu_market;        
sum1 = sharesum_market*numer;
sum11 = 1./(0+sum1);                      
denom1 = sum11(marketForProducts_market,:);    
simShare_true = numer.*denom1;              
EstShare_true = mean(simShare_true,2); 


Y_market = Y_k(ttt,:);

sum_s = simShare_true;
sum_s2 = simShare_true.^2;
alpha_i = alpha_price + alpha_1*Y_market;
sum_alpha_s = alpha_i.*simShare_true;
sum_alpha_s2 = alpha_i.*(simShare_true.^2);

index_d = data_market(:, 3) == data_market(:, 4); 
population_log_d = population_log_k_market(index_d,1);
distance_d = distance_k_market(index_d,1);
FE_d = FE_k_market(index_d,:);
Y_d = Y_market;
xi_d = xi_k_market(index_d,:);
q_0 = q_k_market(index_d,:);

if size(pos_dom,1) > 0

% Create Start Values   
price_d = ones(N_d, 1) * price_k_market(end,1);
q_d = ones(N_d, 1) * q_0;
m1_d = 0.5;
mc_d = 0.9 * price_d;
prod_pct_d = prod_pct(1:N_d,1);

x0 = [price_d; q_d; mc_d; m1_d];
x = x0;

options_fsolve = optimoptions('fsolve', 'Algorithm', 'levenberg-marquardt', 'MaxFunEvals', 50000, 'MaxIter', 5000, 'TolFun', 1e-8);
y_init = solve_mc(x, N_d, alpha_constant, xi_d, alpha_price, alpha_population, population_log_d, alpha_distance, distance_d, q_0, delta_market, index_d, alpha_1, Y_d, FE_d, alpha_FE, nn, expmu_market, Y_market, price_k_market, q_k_market, prod_pct_d);
if any(isnan(y_init)) == 1
while any(isnan(y_init)) == 1 
price_d = rand(1) * ones(N_d, 1) * price_k_market(end,1);
m1_d = rand(1) * 10;
mc_d = rand(1) * price_d;
x0 = [price_d; q_d; mc_d; m1_d];
x = x0;  
y_init = solve_mc(x, N_d, alpha_constant, xi_d, alpha_price, alpha_population, population_log_d, alpha_distance, distance_d, q_0, delta_market, index_d, alpha_1, Y_d, FE_d, alpha_FE, nn, expmu_market, Y_market, price_k_market, q_k_market, prod_pct_d);
end
end
solve_mc_handle = @(x)solve_mc(x, N_d, alpha_constant, xi_d, alpha_price, alpha_population, population_log_d, alpha_distance, distance_d, q_0, delta_market, index_d, alpha_1, Y_d, FE_d, alpha_FE, nn, expmu_market, Y_market, price_k_market, q_k_market, prod_pct_d);  
[x_opt, f_opt, p_flag] = fsolve(solve_mc_handle, x0, options_fsolve);


if (p_flag == 0 || p_flag == -2 || p_flag == 3 || p_flag == 4)
i = 0;
p_flag_new = -2;
while ((p_flag_new == 0 || p_flag_new == -2 || p_flag_new == 3 || p_flag_new == 4) && i < 30)
i = i + 1   
price_d = rand(1) * ones(N_d, 1) * price_k_market(end,1);
m1_d = rand(1) * 10;
mc_d = rand(1) * price_d;
x0 = [price_d; q_d; mc_d; m1_d];
x = x0;
y_init = solve_mc(x, N_d, alpha_constant, xi_d, alpha_price, alpha_population, population_log_d, alpha_distance, distance_d, q_0, delta_market, index_d, alpha_1, Y_d, FE_d, alpha_FE, nn, expmu_market, Y_market, price_k_market, q_k_market, prod_pct_d);
if any(isnan(y_init)) == 1
while any(isnan(y_init)) == 1 
price_d = rand(1) * ones(N_d, 1) * price_k_market(end,1);
m1_d = rand(1) * 10;
mc_d = rand(1) * price_d;
x0 = [price_d; q_d; mc_d; m1_d];
x = x0;  
y_init = solve_mc(x, N_d, alpha_constant, xi_d, alpha_price, alpha_population, population_log_d, alpha_distance, distance_d, q_0, delta_market, index_d, alpha_1, Y_d, FE_d, alpha_FE, nn, expmu_market, Y_market, price_k_market, q_k_market, prod_pct_d);
end
end
solve_mc_handle = @(x)solve_mc(x, N_d, alpha_constant, xi_d, alpha_price, alpha_population, population_log_d, alpha_distance, distance_d, q_0, delta_market, index_d, alpha_1, Y_d, FE_d, alpha_FE, nn, expmu_market, Y_market, price_k_market, q_k_market, prod_pct_d);  
[x_opt_new, f_opt_new, p_flag_new] = fsolve(solve_mc_handle, x0, options_fsolve);
end

if p_flag_new == 1
    x_opt = x_opt_new;
    f_opt = f_opt_new;  
    p_flag = p_flag_new;
end    
end
    
x = x_opt;

price_d = x_opt(1:N_d, 1); 
q_d = x_opt((N_d + 1):(2*N_d), 1); 
mc_d = x_opt((2*N_d + 1):(3*N_d), 1);
m1_d_1 = x_opt(end, 1); 
m1_d = ones(N_d, 1) * m1_d_1;

data_k_ttt = data_market(1, 1:3);
data_k_summary = [data_k_summary; data_k_ttt];

[share_d, profit_d] = get_share(x_opt, N_d, alpha_constant, xi_d, alpha_price, alpha_population, population_log_d, alpha_distance, distance_d, q_0, delta_market, index_d, alpha_1, Y_d, FE_d, alpha_FE, nn, expmu_market, Y_market, price_k_market, q_k_market, prod_pct_d);
f_ttt = profit_d(end, 1);
f_summary =  [f_summary; f_ttt];

price_ttt = price_k_market;
price_ttt = price_ttt(index_d == 0, 1);
price_ttt = [price_ttt; price_d];
price_k_summary = [price_k_summary; price_ttt];

mc_ttt = 0 * price_k_market;
mc_ttt = mc_ttt(index_d == 0, 1);
mc_ttt = [mc_ttt; mc_d];
mc_summary = [mc_summary; mc_ttt];

q_ttt = q_k_market;
q_ttt = q_ttt(index_d == 0, 1);
q_ttt = [q_ttt; q_d];
q_k_summary = [q_k_summary; q_ttt];

share_ttt = share_k_market;
share_ttt = share_ttt(index_d == 0, 1);
share_ttt = [share_ttt; share_d];
share_k_summary = [share_k_summary; share_ttt];

distance_ttt = distance_k_market;
distance_ttt = distance_ttt(index_d == 0, 1);
distance_ttt = [distance_ttt; distance_d * ones(N_d, 1)];
distance_k_summary = [distance_k_summary; distance_ttt];

population_log_ttt = population_log_k_market;
population_log_ttt = population_log_ttt(index_d == 0, 1);
population_log_ttt = [population_log_ttt; population_log_d * ones(N_d, 1)];
population_log_k_summary = [population_log_k_summary; population_log_ttt];

FE_ttt = FE_k_market;
FE_ttt = FE_ttt(index_d == 0, :);
FE_ttt = [FE_ttt; ones(N_d, 1) * FE_d];
FE_k_summary = [FE_k_summary; FE_ttt];

m1_d_ttt = m1_d_1 * ones(size(price_ttt,1), 1);
m1_summary = [m1_summary; m1_d_ttt];

N_d_ttt = N_d * ones(size(price_ttt,1), 1);
N_true_summary = [N_true_summary; N_d_ttt];

market_id_ttt = ttt;
market_id_summary = [market_id_summary; market_id_ttt];

p_flag_ttt = p_flag;
p_flag_summary = [p_flag_summary; p_flag_ttt];

else

data_k_ttt = data_market(1, 1:3);
data_k_summary = [data_k_summary; data_k_ttt];
  
    
%f_ttt = 0 * ones(size(data_k_ttt,1), 1);
f_ttt = 0;
f_summary =  [f_summary; f_ttt];

price_ttt = price_k_market;
price_k_summary = [price_k_summary; price_ttt];    

mc_ttt = 0 * price_k_market;
mc_summary = [mc_summary; mc_ttt];   

q_ttt = q_k_market;
q_k_summary = [q_k_summary; q_ttt];    

share_ttt = share_k_market;
share_k_summary = [share_k_summary; share_ttt];

distance_ttt = distance_k_market;
distance_k_summary = [distance_k_summary; distance_ttt];

population_log_ttt = population_log_k_market;
population_log_k_summary = [population_log_k_summary; population_log_ttt];

FE_ttt = FE_k_market;
FE_k_summary = [FE_k_summary; FE_ttt];

m1_d_ttt = 0 * ones(size(price_ttt,1), 1);
m1_summary = [m1_summary; m1_d_ttt];

N_d_ttt = N_d * ones(size(price_ttt,1), 1);
N_true_summary = [N_true_summary; N_d_ttt];

market_id_ttt = ttt;
market_id_summary = [market_id_summary; market_id_ttt];

p_flag_ttt = 10;
p_flag_summary = [p_flag_summary; p_flag_ttt];

end

clearvars -except i_k ttt data_k_summary f_summary price_k_summary mc_summary q_k_summary share_k_summary distance_k_summary population_log_k_summary FE_k_summary m1_summary N_true_summary market_id_summary p_flag_summary market_id_k data_k prodsMarket companies_k share_k price_k q_k FE_k population_log_k distance_k xi_k alpha_constant alpha_price alpha_population alpha_distance alpha_1 alpha_FE Y_k nn prod_pct

end


filename_f_opt_matlab = sprintf('cost_estimates/data_f_opt_%i_k_test.mat', i_k);  
data_f_opt = array2table([market_id_summary, data_k_summary, f_summary, p_flag_summary]);
data_f_opt.Properties.VariableNames = {'market_id', 'Year', 'Quarter', 'Declarant', 'f_opt', 'p_flag'};  
save(filename_f_opt_matlab, 'data_f_opt');  


